Tight binding models for ultracold atoms in honeycomb optical lattices 
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We discuss how to construct tight-binding models for ultra cold atoms in honeycomb potentials, 
by means of the maximally localized Wannier functions (MLWFs) for composite bands introduced 
by Marzari and Vanderbilt [T]. In particular, we work out the model with up to third-nearest 
neighbors, and provide explicit calculations of the MLWFs and of the tunneling coefficients for the 
graphene-lyke potential with two degenerate minima per unit cell. Finally, we discuss the degree of 
accuracy in reproducing the exact Bloch spectrum of different tight-binding approximations, in a 
range of typical experimental parameters. 



Introduction. Ultracold atoms in optical lattices are 
routinely employed as simulator of condensed matter 
physics, thanks to the possibility of engineering several 
geometry configurations and of tuning the system param- 
eters with great flexibility and precision [21 [3] . In par- 
ticular, honeycomb lattices are attracting an increasing 
interest as they allow to mimic the physics of graphcnc, 
where the presence of topological defects in momentum 
space, the so-called Dirac points, leads to remarkable rel- 
ativistic effects [1HTU]. 

Despite the fact that the potentials describing optical 
lattices are continuous and can be expressed in simple 
analytic forms as the combination of a number of sinu- 
soidal potentials, from the theoretical point of view it 
is often convenient to describe the system by means of 
a tight binding approach on a discrete lattice. In fact, 
the potential intensity can be tuned to sufficiently high 
values so to localize the atoms in the lowest vibrational 
states of the potential wells, justifying a description in 
terms of tunneling coefficients related to the hopping be- 
tween neighboring sites (the potential minima), and in- 
teraction strengths which characterize the onsite interac- 
tion among the atoms [2j. This applies also for the case 
of honeycomb lattices, where a number of tight binding 
approaches have been considered recently [5J [HI EH E] , 
in analogy to the case of graphene [T2Tll5| . 

A crucial ingredient for the connection between the 
continuos and discrete versions of the system hamiltonian 
is the existence of a basis of functions localized around 
the potential minima. This is important not only con- 
ceptually - in order to justify the tight binding expansion 
- but also from the practical point of view, as a precise 
knowledge of the basis functions is needed to connect the 
tight binding coefficients with the actual experimental 
parameters [T5]. In the case of optical lattices with a 
cubic-like arrangement (with a single well per unit cell) 
this basis is provided by the exponentially decaying Wan- 
nier functions discussed by Kohn (2[ [TTl [18] , from which 
one can derive analytic expressions for the tight-binding 



coefficients |19j . However, in general this approach may 
fail when the potential has more than one well per unit 
cell. For example, for the case of a honeycomb potential 
with two degenerate minima in the unit cell, the Kohn- 
Wannier cannot be associated to a single lattice site as 
they occupy both cells for symmetry reasons JT3J H] • 

A powerful approach, that is widely used for describ- 
ing real material structures, is represented by the maxi- 
mally localized Wannier functions (MLWFs) introduced 
by Marzari and Vanderbilt pQ . The MLWFs are obtained 
by minimizing the spread of a set of generalized Wannier 
functions by means of a suitable gauge transformation of 
the Bloch cigcnfunctions for a composite band. Due to 
their exponential decay [H] , the MLWFs provide an 
optimal basis set for tight-binding models and is widely 
employed in condensed matter physics |22) . 

In this paper we discuss the tight-binding expansion 
up to third-nearest neighbors for ultracold atoms in hon- 
eycomb lattices, and explicitly calculate the MLWFs and 
the tunneling coefficients for a potential with two degen- 
erate minima in the unit cell, by using the WANNIER90 
package [23] ■ For the tunneling coefficients we also pro- 
vide an analytic expression in terms of the lattice inten- 
sity, obtained from a fit of the data. Then we discuss 
the validity of different tight-binding approximations - 
including only the nearest neighbor tunneling or up to 
the third-nearest neighbor - in terms of the experimental 
parameters. 

Tight binding approach for the honeycomb poten- 
tial. Let us start by considering the two-dimensional 
graphene-like lattice discussed by Lee et al. [5] 

V(r) = 3 + 2 cos [(&! + b 2 )-r] + 2j2 cos ( b i ' r ) W 

1=1,2 

where r = (x,y) and b 1/2 = V3k L (e x qp v / 3e y )/2 (k L 
being the modulus of the laser wavevectors), correspond- 
ing to a honeycomb structure, with a diamond-shaped 
elementary cell with basis A and B as shown in Fig. [T] 
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FIG. 1. (Color online) Sketch of the honeycomb lattice struc- 
ture and the diamond-shaped elementary cell with basis A and 
B. The length of each side of the hexagon is a = 4n / (3y/3kL) ■ 
The different tunnelings defined in the text are indicated for 
a site of type A. 



The Bravais lattice B in real space is generated by two 
fundamental vectors ai , Si2 defined by a,; ■ hj — 2ir5ij , so 

that B = {jiax + j 2 3L 2 \ji, h = 0- ±1, ±2 ■ • •} 0. 

The starting point for constructing a tight-binding 
model is the many-body hamiltonian for bosonic or 
fermionic particles, described by the field operator 4>{r). 
In the following we will focus on the single-particle term 



dr $ (r)Hoi/j(r) 



(2) 



with Sb indicating the first Brillouin zone, and U S 
U(N) being a gauge transformation that minimizes 
the Marzari-Vanderbilt localization functional £1 — 
Si/ [( r2 )v — ( r )V\ PQ- The real character and the expo- 
nential decay of MLWF functions has been demonstrated 
under very general assumptions [20l EI]- Thus, these 
functions represent an optimal basis for tight-binding 
models. The real character of the calculated Wannier 
functions has been confirmed in our computations. 

In our case, since the honeycomb unit cell contains two 
potential minima, A and B, we can construct a basis of 
MLWFs by considering a composite band consisting of 
the two lowest Bloch bands, that is N = 2. The men- 
tioned Bloch sub-bands have been obtained with a modi- 
fied version of the QUANTUM-ESPRESSO package [25] 
intended to solve the single particle Schrodinger equation 
associated to the external potential |l]). We consider a 
plane wave expansion of the Bloch states, reaching con- 
vergence with an energy cutoff corresponding to £",,=10.5 
En and a k-point mesh of 14 x 14. As a next step, the ML- 
WFs have been computed considering the WANNIER90 
program [22l E3] . The typical shape of the calculated 
MLWFs is shown in Fig. [2] 

The strong localization of |woy(^)| 2 (y — A, B) around 
sites A and B, and their exponential decay are clearly 
visible in panel (a). Panel (c) shows the distribution of 
l w 0A( r )| 2 around the original unit cell j = 0; the figure 
reveals an appreciable overlap of the MLWF with the 



as the mapping onto the tight-binding model is deter- 
mined by the spectrum of the single particle hamilto- 
nian H a = -(h 2 /2m)V 2 + sE R V(r) [21] . Here s repre- 
sents the potential amplitude in units of the recoil energy 
E R = h 2 k 2 L /2m. 

In general, when the potential wells are deep enough, 
the hamiltonian ^ can be conveniently mapped onto a 
tight-binding model defined on the discrete lattice corre- 
sponding to the potential minima, by expanding the field 
operator in terms of a set of functions {wjv(r)} localized 
at each minimum, as 



V>(r) 



(3) 



where j = (ji, J2) labels the cell and v is a band index. 
In Eq. (jij) fity {flj v ) represent the creation (destruction) 
operators of a single particle in the cell j, and satisfy the 



usual commutation rules [aj,, 



it 



(following 



from those for the field tji). 

In order to construct a basis of localized functions at 
each site of the lattice here we consider the MLWFs for a 
composite band discussed by Marzari and Vanderbilt pQ . 
These are a set of generalized Wannier functions Wj u , 
defined from a linear combination of Bloch eigenstates 
V>nfc, namely [HHl] 



1 



fSB 



N 



dk e- lkR > £W*0Vw(r) (4) 




(t>) 


/ \ 











FIG. 2. (Color online) Example of the calculated MLWFs for 
s = 15. (a) Profile of the MWLWs |woA(r)j 2 (solid, blue) and 
\woB(r)\ 2 (dashed, red) along the line joining the A and B 
sites (y = 0) in the original unit cell, (b) Profile of woa(t') 
and Wos(^) along the same path as in (a) with a zoom into 
the small values of the MLWFs. Note that woA(r) (wob(t)) 
becomes negative in the neighborhood of site B (A), (c) Con- 
tour plot of the function log |woA(f)| 2 - The solid and dashed 
lines depict the original unit cell and honeycomb lattice, re- 
spectively. See text for explanations. 
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neighboring B and A sites, indicated respectively by the 
yellow and red arrows in Fig. [I] In addition, the MLWFs 
are characterized by nodes in passing from sites of type 
A to B, corresponding to a change of sign (see panel (b)). 

The approach followed here of including two Bloch 
bands is the minimal approximation, and corresponds 
to the generalization of the single band approximation 
usually employed for cubic- type lattices [2]. Within this 
approximation, the hamiltonian ([2| can be written as 



E 

vv'=A,B 



E a Jv fi 



'3 v 



v 3 „ 



H \wj>„>} 



njtb 



(5) 



where the expansion coefficients in the above equation 
correspond to the onsite energies E v — (wj^H^Wjv), 
and to the tunneling amplitudes between different sites. 
Both quantities are real thanks to the properties of the 
MLWFs. Here we will include up to the third-nearest 
neighbor tunnelings T*^, = — (wj I/ \Ho\w/j^.i-\ l ,/), with 
i = (0, ±1;0,±1) (they depend only on the relative dis- 
tance owing to the uniformity of the lattice). Then, the 
tunnelings can be divided in three classes, as shown in 
Fig. [l]for the case v — A. 

(i) Terms between A(B) and the three nearest neigh- 
bors of type B(A), yellow arrows in Fig. [I] e.g. 



u jA \H \w jB ) 



(6) 



(ii) Terms t„ between sites of the same type {A or B) 
within neighboring cells, red arrows in Fig. [I] e.g. 



( w {3\+-i-,3*)i'\ H o\ w {ju3z>) 



(7) 



In general the two tunneling coefficients t A and t B are 
different. As a specific example, here we will explicitly 
compute them for the case of degenerate minima where 
t A = t B . In this case we can set t± = — t v , where the sign 
is chosen in order to have t\ positive defined (see Fig. [2] 
and the discussion about the sign of the MLWFs). 

(Hi) Terms connecting A(B) to B(A) at opposite cor- 
ners of the hexagon, blue arrows in Fig. [lj e.g. 



( w UiJ2)a\H \w Ui . 



■1,J2-1)B/ 



(8) 



We remark that the above derivation of the tight- 
binding model is valid in general for any potential with a 
honeycomb structure, with two minima per unit cell, and 
not just for the potential with two degenerate minima in 
Eq. ([I]). In the following we will consider explicitly the 
latter case in order to provide a specific example. 

The behavior of the different tunneling coefficients as 
a function of the lattice intensity s is shown in Fig. [3] 
In order to extract from the numerical values an analytic 
formula we consider a fit of the type tj = As a e~^^ 
(i = 0,1,2), in the range s > 3, with A, a, and j3 as 
fitting parameters. For to we find 



tn 



1.16s 



0.95 -1.634/1 



(9) 



that has to be compared with the semiclassical estimate 
of Lee et al. [5], |t | = 1.861s°- 75 e- 1 - 582 ^ (see Eq. (38) 




FIG. 3. (Color online) Behavior of the various tunnelings as a 
function of the lattice intensity s. The lines are the result of a 
fit of the numerical data, and those for to and ti coincide with 
that extracted from a fit of the Bloch spectrum (see text). 



in [S]). In the range of s considered here, we find that 
the latter overestimates the actual value in ([9| by about 
8% for s = 30, up to about 40% for s = 3. For the other 
two terms we get 



h = 0.78s 1 - 85 e - 3 - 404 ^, 
t 2 = 1.81s 2 - 75 e- 5 - 196 ^. 



(10) 
(11) 



These three fit arc shown as dashed lines in Fig. [3] 

Tight-binding spectrum. A convenient way to check the 
regime of validity of a given tight binding approximation 
is to compare its prediction for the energy spectrum with 
the exact Bloch spectrum. The latter can be readily com- 
puted by means of a standard Fourier decomposition [S] . 
The typical structure of the two lowest bands E± (fe) is 
shown in Fig. [4j It is characterized by Dirac points at 
the vertices of the Brillouin zone (a regular hexagon), 
where the local dispersion is linear and the two bands 
are degenerate [5|. For convenience, we fix E±(ko) = 
at the Dirac points. 

The tight-binding spectrum can be derived as follows 
[IE]. By defining b uk = (^/27r) ^ e 



~ %k ' R ihj v the 



hamiltonian in Eq. 



5h can be written as 



lb 



H 



with h vv i(k) 



E 

j'=A,B 



S,3 



dk h vv ,{kjiy vk b v ik 



(12) 



h V i,i (k) turns out of the form [S] 



^2 i e lkRi T^ v ,. Finally, the matrix 



h vv >{k) 



e A (k) z(k) 
z*(k) e B (k) 



(13) 



As for the diagonal terms, the leading contribution 
comes from the onsite energies E u , that can be conve- 
niently written as E A /b = i £ > by shifting the total en- 
ergy by an overall constant. In addition, there is a cor- 
rection due to t„, so that eventually we get 



£A/s(fc) = E A/B - t A/B F(k), 



(14) 
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FIG. 4. (Color online) (a) Bloch spectrum E±(k) of the low- 
est two bands for s = 5. The hexagon represents the first 
Brillouin zone. (b,c) Cut for k y = and k x = 0, respectively 
(blue, solid lines). The latter are compared with the predic- 
tion of the full tight-binding model (red, dashed line), that 
with to and t\ (black, dotted), and with just to (magenta, 
dot-dashed). Note the asymmetry of the two bands. 

F(k) = 2 cos (k ■ (en - a 2 )) + 2 cos (fe ■ a,) . (15) 

Instead, the off-diagonal terms get the leading contri- 
bution from the term proportional to t [5], and a cor- 
rection due to t 2 , z(k) = t Z (k) + t 2 Z 2 (k) 7 with 

Z (k) = 1 + e- lkai + e- lka2 

Z 2 {k) = e ~ ife '( a l+ a 2) -|- g-*fe'(«l-"2) _|_ g-ife-(o 2 -Ol) 

Finally, by diagonalizing the matrix h vv i (fe) and defin- 
ing t± = (i,A ± ts)/2, we get the following expression for 
the tight-binding spectrum 

e± (fe) - -t+F{k) ± yJ(t_F(k)-e) 2 + \z(k)\ 2 . (16) 



t +h = (e+(0)-e_(0))/6andii = (e+(0) + e_(0))/18. 
In addition, since t 2 is negligible with respect to to ( see 
Fig. [3]), practically the former expression can be used as 
an estimate of to- Notably, these estimates coincide with 
the dashed lines shown in Fig. [3J therefore providing an 
independent check of the tunnelings calculated by means 
of the MLWFs. 

In general, one can evaluate the degree of accuracy 
in reproducing the exact Bloch spectrum by defining an 
energy mismatch as follows 

Se» = ~ dk[E n (k)-e n (k)} 2 (18) 

with AE n being the n— th bandwidth (n — 1,2). The 
results are shown in Fig. [5j This figure shows that the 
tight-binding model with up to third-nearest neighbors 
accurately reproduces the band structure for s > 3, with 
an error below 1%. In fact, this is a range of lattice in- 
tensities where one would expect the MLWFs to localize 
strongly around each minimum (that is, a proper tight- 
binding regime) . Then, while the inclusion of t 2 provides 
only a minor correction, the model with just the nearest 
neighbor tunneling to [5] is clearly less accurate, reaching 
the level 5e n < 1% only for s > 15. This may be par- 
ticularly relevant for the range of parameters of current 
experiments [S]. 
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Again, this is valid for a generic honeycomb structure 
with two minima per unit cell. For the particular case 
of degenerate minima, as for the potential in 0, this 
expression further reduces to 

e±(fe) = hF(k) ± \t Z (k) + t 2 Z 2 {k)\ + 3ti (17) 



where the last term has been added in order to make 
the energy vanishing at the Dirac points, consistently 
with the definition of the Bloch spectrum. A specific 
example is shown in Fig. |4^b,c), for s = 5. The figure 
shows that the tight-binding model with just to is not 
sufficient to reproduce the band structure, and that at 
least the inclusion of t\ is needed. In particular, the latter 
is necessary to account for the band asymmetry, see Eq. 
([T7|). It also remarkable to notice that the values of t a +t 2 



and t\ can be extracted from a fit of the Bloch spectrum 
at k = 0, by using the expression (17). In fact, we have 



FIG. 5. (Color online) Energy mismatch Se n for the first (a) 
and second (b) band, for different levels of approximation of 
the tight-binding model. 

Conclusions. In this article we have demonstrated the 
power of the maximally localized Wannier functions for 
composite bands in determining the parameters of tight- 
binding hamiltonians describing ultracold atoms in opti- 
cal lattices. The application to a honeycomb structure, 
directly connected to the graphene physics, allows us to 
accurately parametrize the optimal tight-binding param- 
eters, providing a thorough analysis of the range of va- 
lidity of different approaches. 
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